Multiple barriers in forced rupture of protein complexes 
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Curvatures in the most probable rupture force (/*) versus log-loading rate (logr/) observed in dy- 
namic force spectroscopy (DFS) on biomolecular complexes are interpreted using a one-dimensional 
free energy profile with multiple barriers or a single barrier with force-dependent transition state. 
Here, we provide a criterion to select one scenario over another. If the rupture dynamics occurs 
by crossing a single barrier in a physical free energy profile describing unbinding, the exponent i/, 
from (1 — f* I fc)^^" ~ (logr/) with /c being a critical force in the absence of force, is restricted to 
0.5 < < 1. For biotin-ligand complexes and leukocyte-associated antigen-1 bound to intercellular 
adhesion molecules, which display large curvature in the DFS data, fits to experimental data yield 
V < 0.5, suggesting that ligand unbinding is associated with multiple-barrier crossing. 
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INTRODUCTION 

Single molecule pulling experiments have generated a 
wealth of data, which can be used to probe aspects of 
folding that were not previously possible [TH3]. In ad- 
dition, DFS has been used to decipher the energy land- 
scape of molecular complexes by measuring the rupture 
force (/) by linearly increasing load at a rate rf (= 
df/dt). Because of the stochastic nature of the unbind- 
ing events, / varies from one complex (or realization) 
to another, giving rise to an r/-dependent rupture force 
distribution (P(/)). For a molecular complex obeying 
Bell's formula, k{f) — koffexp{fx^/kBT), Evans and 
Ritchie showed that the most probable force is /* — 
{kBT/x^)\og{rfX^ /koffkBT) g], where x^{= Xts-Xb) is 
the location of the transition state (xts) from the bound 
state (xf,) projected along the pulhng coordinate and koff 
is the unbinding rate in the absence of force. However, 
Bell's formula is applicable only if the molecular com- 
plexes are mechanically brittle or if the applied tension 
is sufficiently small that x^ does not shift upon applica- 
tion of force [5]. More generally, /* follows a (logr/)"^ 
dependence where v depends on the details of the as- 
sumed one dimensional (ID) model potential [5til3j. The 
basic assumption in all these works is that a single free 
energy barrier along the pulling coordinate is sufficient 
to describe force-driven rupture of the bound complex. 

Sometime ago Merkel et al. used DFS to probe load de- 
pendent strength of biotin bound to ligands, streptavidin 
and avidin [14] , showing that over six orders of variation 
in rf (from about 10^^ to rf in excess of 10^ pN/s) the 
plot of /* versus logr/ ([/*,logr/] plot) varies nonlin- 
early for both ligands. We note parenthetically that it 
is also common to observe curvature in unfolding rates 
of proteins when the r/ is varied |15j . By careful data 
analysis combined with molecular dynamics simulations 



they proposed an energy landscape for the complex, with 
multiple energy barriers jl4j. A similar picture emerges 
in the rupture of intercellular adhesion molecules (ICAM- 
1 and ICAM-2) bound to leukocyte function-associated 
antigen-1 (LFA-1) upon application of force (THj. 

In principle, however, nonlinearity in [/*,logr/] plot 
could also arise from load dependent variation in x^ 
[17] in a ID energy landscape with a single barrier [5]- 
[TOl fT2| [T3l [TT] . a theoretical model describing force- 
induced escape from a bound state with a single barrier 
in a cubic potential {v = 2/3) has been used to ratio- 
nalize the biotin-ligand data by identifying various linear 
regimes demarcated by r/ 0. However, in the absence of 
easily discernible changes in the slopes in [/*,logr/] plot 
it is difhcult to justify such an analysis. Here, we show 
by analyzing experimental data that the observed non- 
linearity in the DFS data of several protein complexes can 
be better accounted for with an energy landscape con- 
taining multiple sequential barriers, as originally demon- 
strated [li l fB ] . 

To illustrate how steep curvatures in DFS data can 
arise naturally from a ID free energy profile we cal- 
culated P{f) and [/*,logr/] of forced-escape kinetics 
of a quasiparticle from a potential with two barriers, 
U{x) = Ax{x - l)[{x -2){x- 3)(x - 4)(x - 5) -t- 1] with 
A > (Fig[T]). The distributions P{f) are typical of 
what is observed in experiments (Fig[I]3). For all val- 
ues of A, [/*, logrj] plots are curved although one could 
discern a modest change in slope (FigjlJ:). The load- 
ing rate dependent x^{rf)^ calculated from the slope of 
the data ksT /x^{r f) at each rf in FigjlJ;, changes from 
^ 3 nm to < 1 nm. The precipitous change in x^ at 
r/ K, [e~^ — e°) pN/s reflects the transition of the confin- 
ing barrier from outer to inner barrier with an increasing 
force (see FigjT^). In contrast, gradual change of x^ in 
the range < logr/ < 10 is most likely due to the move- 
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FIG. 1. Rupture characteristics obtained numerically using 
a potential with two barriers at constant loading rates, (a) 
[/(a;) (magenta) and Uix) — f ■ x (cyan) with A — 5 pN ■ nm 
and / = 50 pN. Reflecting and absorbing boundary condi- 
tions are set at a; = a and x — b, respectively, (b) Rupture 

force distributions, P{f) = k{f)/rf ■ exp [- /(f df'fc(/')/y/] , 
at varying r/ were computed by using mean first passage time 
(MFPT), k-\f) = D-i/;dye«^(«)--f-«)/Jd2e-«^(^)-^-^), 
starting from the first bound state at a{—0 nm) to reach an 
absorbing boundary at 6(=5 nm). MFPT expression is valid 
in the force regime where stationary flux approximation holds 



[4]. The length was scaled by nm, and D = 1.0 x 10^ 



7^ 



was used for the diffusion constant, (c) [/*,logr/] plots at 
three A values. Fits of [/*,logr/] to Eqis] yield u <C 0.5 for 
all A values {u = 0.064, 0.075, 0.046 for A = 4, 5, 6 piV ■ nm, 
respectively). In this case, the data should be divided into 
two regions and analyzed by the two linear fits as depicted 
using green lines on the curve with A = 6pN ■ nm. (d) load- 
ing rate dependent x^{rf){= xts ~ Xb), extracted from the 
slope of plot at each r/ in (c) with A — 5 pN ■ nm, shows a 
sharp decrease from ~3 nm to < 1 nm around r/ ~ (e^^ "g") 
pN/s. 



ment of the inner transition state (see Fig.5C in Ref.(5)). 
Although it is straightforward to interpret that the two 
discrete slopes in Figjlj: (or the precipitous transition of 
in Fig[l|l) are due to crossing two barriers since the 
underlying potential is given in Figjl^, it is nontrivial to 
solve the inverse problem of unambiguously determining 
from [/*,logrj] plots and decide whether the underly- 
ing free energy profile has a single barrier with a moving 
transition state as increases or multiple barriers. 

To establish a criterion for ascertaining whether the 
energy landscapes for forced-ligand rupture from biotin 
and LFA-1 have multiple barriers, we study the range of 
applicability of DFS formalism based on a model poten- 
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FIG. 2. The n-dependent shape of G{x) (Eq[T]). The poten- 
tial with increasing n is associated with more brittle molecular 
complexes. The yellow circle {x = Xc) denotes an inflection 
point and a cusp in each even and odd n potential, respec- 
tively. 



tial with a single barrier. Consider a Kramers' problem 
of barrier crossing in a free energy profile G{x) in which 
a single barrier separates the bound and unbound states 
of a quasi-particle as in ligand bound in a pocket of a 
receptor: 

G{x) = G{x,) + Ux- Xe) + \ [x - x,)"+i (1) 

[n + I)! 

with M > 0. In G{x), a ID free energy profile with a sin- 
gle barrier, the shape of barrier and energy well is approx- 
imated using n-th order polynomial with n = 1, 2, 3, • • • . 
For odd n, we assume that G{x) — ~oo for x > Xc, so 
that the transition state of G{x) is cusped. In the ab- 
sence of tension, the barrier height, G^, and the location 
of transition state, x^, are G^ = X^^^fcin-lfc/My/" and 

= xi''^Uc/My^^ : respectively, where x — 1 (foi' odd 
n), 2 (for even n). Thus f, = '^G^/x^. The form 
of G{x), an extension of the microscopic theories using 
harmonic-cusp or linear-cubic potential, accounts for the 
degree of plasticity (or ductility) or brittlcness of the 
energy landscape by changing n (Fig(2]) [5.. Under 
tension, Geff{x) — G{x) ~ f ■ x; fc should be replaced 
with /,(! - ///,) = U. Therefore, GHf) = G^ei+i/" 
and x^if) = xtei/". Although Eq{l] looks similar to 
the one Lin et al. used to discuss rupture dynamics 
for £ ^ 1 where the barrier height is almost negligi- 
ble |12j . we did not impose any specific force condition 
on G{x). Instead of attributing the movement of tran- 
sition state to a large external tension [THTUl [T^ [T^ . 
we mapped the non-linearity in DFS data onto G{x) 
that has the n-dependent shape of transition barrier and 
bound state. In G{x), increasing brittleness makes x^{f) 
insensitive to applied tension, which is dictated by n; 



x^{f)/x^ = e^/" 1. For a generic free energy pro- 
file F(x) with high curvatures at both x = Xts and Xf,, 
xi{f)/xi = 1 - J/x^ X {\F"{xts)\-' + \F"{x,)\-^) ^ 1 
5J. When free energy profile is associated with a brittle 
barrier, Bell's formula can be used to extract the feature 
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of the underlying ID profile from DFS data [5]. 

For general n, the Kramers rate equation based on the 
EqjT] under tension can be derived as: 



k{e) 



exp 



(2) 



where k is the prefactor in Kramers theory and a{n) = 
— 1/n) with X = 1, 2 for odd and even n, respectively. 
For a given fc(/), the most probable unbinding force is 
determined by dP{f)/df\f=f* = 0, resulting in a general 
equation for /*: 



which leads to 

^ -1 

log 
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Under the typical condition that rupture occurs by ther- 
mal activation, i.e., / ^ /c(e « 1) and PG^ ^ 1, the 
most probable unbinding force is approximated as: 



1- - 



(5) 



where v — In deriving Eqjsl using G(x), the large 

force e{= 1 — ///c) <C 1 or fast loading condition, an 
assumption made in obtaining the mean unbinding force 
expression similar to Eqjs] [6l [T2j [13] , is not needed. Only 
the shape of the energy potential matters in deriving EqjS] 
from EqjT] The DFS data will have a larger curvature for 
smaller n, namely when the energy landscape associated 
with a protein complex is more ductile (Fig. 2). Because 
n = \ (harmonic cusp), 2 (linear cubic), . . ., cx) (Bell), v 
must satisfy the bound. 



1/2 < 1/ < 1 



(6) 



for an arbitrary ID profile that suffices to describe rup- 
ture kinetics. 

For forced-rupture of biotin-ligand complex, fits to the 
entire range of the data using Eq|5] give v in the disal- 
lowed range; v — 0.40 (biotin-streptavidin) and v — 0.070 
(biotin-avidin) (see Fig[3^). Even in biotin-streptavidin 
case, the parameters extracted from the fits with v — 0.40 
and V — 0.5 (fixed) are comparable; however, the fit with 
V = 0.40 is superior yielding both smaller relative error 
and reduced chi-square value, Xredi than with v — 0.5, 
the lower bound of Eq. [6] that gives the maximal curva- 
ture in the single-barrier picture (see Figj3]ja) and Fig S2 
in the SI). For both biotin-ligand complexes, our crite- 
rion consistently suggests that the unbinding landscapes 
for the complexes involve more than one barrier, as was 
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FIG. 3. Analysis of DFS data with large curvatures, (a) The 
data obtained using biomembrane force probe (BFP) with 
force constant in the range 0.1-3 pN/nm [T^ were fitted to 
Eqjs] (solid lines) with v=QAQ for biotin-streptavidin (circle) 
and i^=0.070 for biotin-avidin (triangle). The x^[rf) at each 
r/ is calculated on the right using the slope of four succes- 
sive data points of [/*,logr/] plot. Analyses of data using 
restricted v values {v = 0.5 fit is in dashed line in Figjs^) are 
in the SI (b) Analysis of DFS data of LFA-1 and its ligands, 
ICAM-1 and ICAM-2 in Ref. [K]. The fits in log-log scale 
are shown on the right. In all cases, p < 1/2 suggests that 
for these complexes as well the underlying free energy profiles 
must contain at least two barriers; thus multi-state fits are 
required by dividing the DFS data into multiple regions as 
was already surmised in [TB] . 



emphasized by Merkel et al. [H]. Next, we analyzed 
the extensive data on LFA-1 expressed in Jurkat T cells 
whose binding affinity to ICAM-1 and ICAM-2 can be 
enhanced by treating the cells with phorbol myristate ac- 
etate (PMA) and the divalent counterion, Mg^+. Under 
all conditions the exponents that best fit the DFS data 
are i' < 0.5 (FigJSc)). As originally argued by entirely 
different method [ig rupture of ICAM-1 and ICAM-2 
from LFA-1 is best described using a free energy profile 
with at least two barriers. Taken together we arrive at a 
consistent conclusion that v < 0.5 implies that the under- 
lying free energy landscape in protein-ligand complexes 
has multiple barriers. 
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Mathematically the inequality (Eq|6| is not strict be- 
cause it is possible to construct ID profiles with v < 0.5 
for which [/* , log r f] plots exhibit curvature like those ob- 
served in experiments. However, such free energy profiles 
are physically pathological with non-existing first deriva- 
tives in the vicinity of the bound complex and fits to the 
data give manifestly unrealistic parameters (see Support- 
ing Information for detailed calculations and analysis). 
For the physical free energy profiles Eq|6] is rigorously 
satisfied. In addition, there is no compelling reason to 
choose a special n value even if ID profile is deemed ad- 
equate, and thus v ought to be treated as a parameter. 
Although Ref. [TT] used as a free parameter, the va- 
lidity range for i/ was not discussed. If a global fit of 
[/*,logr/] data using Eqj5] yields v < 0.5 and the effect 
of probe stiffness [THj is absent in the DFS data (see be- 
low) , we can conclude that a single barrier description of 
the energy landscape is inadequate. 

In principle curvature in the DFS data could also arise 
due to probe stiffness. Simple procedure of tiling free en- 
ergy profile by the amount —f-x under tension is widely 
used in analyzing single molecule force experiment. How- 
ever, more rigorous formulation for the effective free en- 
ergy profile under load using a transducer with stiffness 
k should read Gtot{x,Xtr) = G{x) + ^kf,ff{x — XtrY 
where x is the position of the end of molecule, Xtr 
is the position of transducer, and fee// is the effective 
stiffness of molecular construct combining the trans- 
ducer and the complex (fcj/y = k^^^ + fc~^). In fact, 
\keff{x - XtrY = -f ■ X + \kf,ffx'^ + ^keffXf^ with 
/ — kef f Xtr- Therefore, the effective free energy for 
the complex under tension should be written in general 
as Geff{x) = G{x) — [f — ^keffx] ■ X [in]. As long as 
/ ^ ^keffX (or Xtr ^ 2;/2) especially when k^ff is 
small as in optical tweezers or BFP, one can approx- 
imate Geff{x) ~ G{x) — f-x. Otherwise, rebinding 
from transient capture well created by a large probe stiff- 
ness at near-equilibrium loading condition could give rise 
to the nonlinearity in the DFS data [T^. Therefore, 
the rupture force being measured should be replaced by 
f* f* - IkeffX^e'^/"-, and at low forces (/ <C fc) the 
most probable force measured by using a transducer with 
high stiffness such as AFM could be approximated as. 
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The effect of probe stiffness manifests itself as a non- 
vanishing plateau force, which could be as large as fpi « 
(10 - 100) pN when k^ff ~ 100 pN/nm and = 0.1 - 1 
nm, even when is small enough that /dfs — [18 . 

The biotin-ligand complexes data preclude this possi- 
bility because the probe stiffness of BFP k^ff = 0.01 — 0.3 



pN/nm |14| . which is 1-2 orders of magnitude smaller 
than the fee/ / value discussed in the literature [201 HI] ■ 
The value of fee// is 0.5 — 2.0 pN/nm in the experiments 
involving LFA-1 [16,. Even the largest estimated x^ value 
for the outmost barrier {x^ « 3 nm) only yields fpi < 1 
pN. Furthermore, if the nonlinear curvature of DFS data 
is suspected to be due to the stiffness effect, this ought 
to be discerned from the curvature due to multiple barri- 
ers by producing DFS data at a reduced probe stiffness. 
The curvature due to multiple barrier should persist in 
the DFS data even with a small fce//- Thus, the cur- 
vature in the data in [14] can only be attributed to the 
presence of multiple barriers. 

The condition (Eqj6| for single-barrier based ID theo- 
ries for DFS |7lfT5] provides a guideline to judge whether 
the curvature in DFS data is due to multiple barriers or 
single barrier with a ductile transition state. Our work, 
which does not consider complications due to various 
multidimensional landscape scenarios ^21122, shows that 
the extracted parameters from the data for the protein 
complexes with ligands using ID profile with multiple 
barriers are physically reasonable [HI [16]. Additional 
justification for the use of such energy landscapes can 
only be made by studying the structures of the protein 
complexes in detail. 
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SUPPORTING INFORMATION 

DPS theory for a free energy profile with non- 
integer n. It could be argued that the inequality 
1/2 < < 1 (valid rigorously for integer n) that accounts 
for the curvature of DFS data is not mathematically re- 
quired. Here we show that it is possible to construct ID 
free energy profiles for which i' is clearly less than 0.5. 
Indeed, these free energy profiles can even have nearly 
vanishing v. However, such profiles are unphysical be- 
cause near the bound state they have incorrect curva- 
tures compared to the physical profiles discussed in the 
text and in the references cited therein. More impor- 
tantly, the first derivatives of these free energy profiles, 
which yield i' < 0.5 do not exist near the bound state 
i.e, they have a singularity. For these and other reasons 
(see below) we reject these free energy profiles as plausi- 
ble models for explaining the curvatures in the observed 
[/*,logr/] plots in a number of protein complexes dis- 
cussed in the main text, which have all been explained 
using a two barrier model. 

A free energy profile with < n < 1 (see Eqjl]) can 
lead to < < 1/2 since ly — To explore this 

possibility, we consider a free energy profile, 



G(x) 



c^/'" - bx 



(SI) 



with a > 1. Here we may regard n = l/tr, and hence 
n < 1. The term —bx is required for the potential to 
have a barrier at a finite value of x^ , the location of the 
transition state (TS). The TS location and the associated 
barrier height are 



...(^)* 

= {<T-l)a 



(S2) 



Under tension / the potential becomes Gefj ~ G{x) 
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and the barrier height are 
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where ry = (1 + //Aat) with fi/„ = Therefore, 
one can rewrite Eq SI {G{x)) and an effective free energy 



{Geff{x)) under tension as 
crGi 
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Given Geff{x), it is possible to obtain the mean first 
passage time expression corresponding to the lifetime of 
the complex that can be measured in single molecule ex- 
periments. It is given by, 

-1 />oo /^X 

k{f)-^ = - da;e^G="(^) / dye-^^^"'^y'> . (S5) 
D Jo Jo 

The saddle-point approximation, k{f) w 

D ^%g?^e-^«*(/) with 

Eq.A3, yields Kramers' equation for the escape rate: 
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with a((7) 



Here, note that the k, de- 



fined as the prefactor, contains the singular integral 

{ibound'^y'^'^^'""'^^^^) ■ ^^^"^g relationship 

= [^(/*)]^/^/ to derive the most probable force, 

we get 
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By assuming fSG'^ ^ 1 and / ^ fi/a, we can simplify 
the above equation into 
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Therefore, the most probable force (/*) for the fractional 
potential (Eq SI ) is 



/* ~ fl/a 
= fl/a 
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where a = was employed in the last line. 

There are a few important comments about Eq|S9| that 
need to be made, (i) Note that the form of Eq |S9| is very 
different from Eqj5j The difference is related to the afore- 
mentioned difficulties associated with Geff{x). Never- 
theless, it is possible mathematically to construct model 
free energy profiles, without regard to physical consid- 
erations, for which [/*,logr/] plot has curvature that is 
reminiscent of what is observed in experiments, (ii) Al- 



though Eq S9 can be used to obtain excellent fits to DFS 
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(logr,.-36.3) 



(b) 



G' = 49.2k„T, .r = 7.51 nm, k-=3.0x10"s ', CT=4.17 
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log r [pN/s] 



/' = 422.1 



1 



-112.6) 



19.1 ' 

G* = 19.1k„T, jc* = 1.42nm, /f=2.68 x IQ-^s ' , (T=l .13 



/• = 11.15 



1 

3n 



(log,;,, -27.5) -1 



G*=31.1k„T, x* = 2.70nm, k:=5.59 X 10"s ' , a=5.29 



FIG. SI. (a) The blue curve is the bare (/ = 0) free energy profile of the form given in Eq. (SI) and the green curve is 
the tilted form of G{x) in the presence of force. By fitting the numerically computed (black circles) /* as a function of r/ to 
Eqjs] (red curve) we obtain the parameters shown below. Although the features of original potential G{x) — lOOa;^''"^ — 10a; are 
reasonably recovered (a is larger than the value in G{x)) by using Eqjsj the extracted value of k is unrealistically large, (b) 
Eqjsjwas used to fit the DFS data of biotin-strepavidin (circles) and biotin-avidin (triangles). Although the quality of fit is 
excellent, the unrealistically large value of k, due to the singularity of the hypothesized fractional potential at a; = 0, suggests 
that the potential with a fractional power should not be used for the analysis. In fact the k values are comparable to or much 
greater than the TST estimate ksT/h (« 6.2 x 10^'^ s"^), which of course makes no sense. Hence, we can rule out the free 
energy profiles of the form given in Eq.Sl to analyze DFS data on protein complexes. 



data on protein complexes, it turns out that the extracted 
value of K is extremely large and are clearly unphysical 
(see Fig SI ). The unphysical values are a consequence of 
the singularity of G{x) (or Geff{x)) at a; = 0. We con- 
clude that the free energy profiles with a fractional power 



of n, which mathematically creates singularity at bound 
state, is not suitable to be used to analyze experimental 
data. Thus, the bound 1/2 < < 1 in Eqj6]must hold 
for physical ID free energy profiles with a single barrier. 
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FIG. S2. Analysis of DFS data using Eq|5] for (a) biotin-streptavidin and (b) biotin-avidin. For each v, the fits, residuals 
(I //it ~ fexp\/flxp X 100), and extracted parameters were summarized in the table on the right. We can draw some general 
conclusions from the fits. For the biotin-streptavidin complex, fit with v = 0.397 produces the smallest errors although at high 
loading rates the relative errors for v — 0.397 and v — 0.5 are comparable. There are variations in other parameters {x^ , G*, 
and k) for all u. For the biotin-avidin complex the situation is far worse. In particular, the relative errors in the fits are large 
even when u is varied. Similarly, the parameters extracted from the fits are not totally consistent. Taken together, the fits 
using a ID free energy profile with a single barrier is not appropriate to describe the rupture kinetics of these two complexes. 



